clear all
close all

% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the limited commitment village economy for
% both the standard specification (individual deviations) and the
% alternative specification (coalitional deviations)

% This version to be used for Shirapur with coalitional deviations

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version September 2021
% =====================================================================


filetosave='some_name'; % name of workspace saved

%====================================================
% select parameters
% ===================================================

% ======== a. individual's income process ========
one_per_aut=0;                     % set to 1 for one period of autarky after deviation; to 0 for immediate risk-sharing
village=3;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=1;                         % set to 1 for coalitional deviations
unconditional=0;                         % only applies to village 1: set this to 0 to load datasheet without unconditional correction, 1 with larges observations eliminated, 2 with first period eliminated altogether
fixedeffect=1;                     % set this to 1 for an income process with fixed effect
ytaxrate=0;                         % linear tax on income
update=0.2;                         %update policy function by this amount
convlength=10;                       % the conv criterion must be fulfilled at least in convlength consecutive periods to avoid jumping 'in'
Tessa1Tobi2=1;
n_rhogrid=25;
estrho1=0;
Nphi=1;
server=0;				% set to 1 to run on server; o wise path for workstations is used
buildup=0;              %set to 1 and change Qrho and Qbeta loop to calculate for all group sizes
waut_restrict=0;             % set to 1 for more restrictive specification of coalitional deviations: individual can offer coaldev only when all but 1 individual in the ROV have high income
Rov_sb=0;
% other execution options
Nincproc=1;                         % number of income processes to look at
T=6200; Tinit=201;                 % number of simulation periods in total, and number of initial periods to discard
maxcerr=0.9;                          % max cons meas error in log levels
N_cerr=30;                           % number of support points for cons meas error
if estrho1==0
    maxincerr=2/3;                     % maximum income measurment error in multiples of income variance (in levels)
elseif estrho1==1
    maxincerr=0;
end
if village==1
    vss_high_inddev=34-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=34-1;
elseif village==2;
    vss_high_inddev=37-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=37-1;
elseif village==3
    vss_high_inddev=31-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=31-1;
end
vssbuildupvec_inddev=[1:9 vss_high_inddev]; % if you choose buildup==1, rhis is the vector of village sizes for which the solution is calculated in the individual deviations case

if Tessa1Tobi2==1
currentfolder='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication';    
elseif Tessa1Tobi2==2
currentfolder='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication';
end
if server==1
    
else
    if Tessa1Tobi2==2
        outputpath=[currentfolder '\Programmes/Main Results and Robustness']
        programpath=[currentfolder '\Programmes']
    elseif Tessa1Tobi2==1
        outputpath=[currentfolder '\Programmes/Main Results and Robustness']
        programpath=[currentfolder '\Programmes']
    end
end

betavec0=[0.5:0.02:0.98];
betavec1=[0.84:0.02:0.88,0.89:0.01:0.99];
betavec3=[0.9:0.01:0.98];
if coaldev<2 || coaldev==3
    rhovec=[1];
elseif coaldev==2
    rhovec=ones(1,Nphi);
end

%main execution options
gapcritset=0.0001; %convergence criterion for change in values


% consumption
rov_average=1;
% set this to 1 if you want to calculate the moments for levels of
% consumption based on log consumption
momentclevinlog=1;

% convergence criterion and maximum number of iterations
maxitgapset=400;
% this is a punishment factor: default involves a consumption penalty in the first period equal
% to caut*Pun_fac
Pun_fac=0.0;
% grid of time-varying planner weights
lambda=[0.05:0.001:0.95]';
m=length(lambda);
if floor((m+1)/2)~=((m+1)/2)
    disp('care: lambda has to have uneven no of points in first dimension')
    stop
end
S_KEEP=nan(1000*3,2,vss_high_inddev);                  % pre-allocate the matrices for coalitional deviations, so that you don't run simulations vss times
S_rov_KEEP=nan(100,1,vss_high_inddev);              % this one is for a maximum of vss=19 for the CD
P_rov_KEEP=nan(100,1000,vss_high_inddev);
PHI_KEEP=nan(m,1000*3,2,vss_high_inddev);
IDIO_DIST_ROV_KEEP=nan(1000*3,3,vss_high_inddev);

% ====================================================
% I. Specify Income process for individual and Rest of Village (rov)
% ====================================================
% data moments for log income from ICRISAT data - all villages
% first col: Aurepalle, secdon: Kanzara, third: Shirapur
if fixedeffect==1
    if unconditional==0
        datafile='armoments_fe0';
    elseif unconditional==1
        datafile='armoments_fe1';
    end
else
    datafile='armoments';
end
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);

if maxincerr>0
    if Nincproc>1
        varnu_max=maxincerr*vardy/2;                % nu is measurement error in log levels
    else
        varnu_max=zeros(1,3);
    end
    
    rho_min=covyylag(village)/vary(village);                      % AR(1) coefficient w/o meas error
    rho_max=covyylag(village)/(vary(village)-varnu_max(village));          % AR(1) coefficient w max meas error
    rho1vec=linspace(rho_min,rho_max,Nincproc);
    varnu=max(vary(village)-covyylag(village)./rho1vec,0); % max prevents the first entry to be -e17
    vareps=(vary(village)-varnu).*(1-rho1vec.^2);
else
    rho1_data=covyylag(village)/vary(village);
    vareps_data=vary(village).*(1-rho1_data.^2);
    rho1vec=[-0.6,-0.4,-0.2,0,0.2,0.4,0.6];
    vareps=vary(village).*(1-rho1vec.^2);
    varnu=zeros(1,length(vareps));                % nu is measurement error in log levels
    Nincproc=length(vareps);
end
% ===================================
% iterate over income processes
% ===================================
disp('now solving for village,coaldev, unconditional, Nincproc, fixedeffect one_per_aut Rov_sb - press any key to continue')
disp([village coaldev unconditional Nincproc fixedeffect one_per_aut Rov_sb])
pause
indicator1=zeros(1,size(rhovec,2));
for incproc=1:Nincproc
    disp(incproc)
    
    if buildup==0 && maxincerr>0 && Rov_sb==0 && one_per_aut==1
        filetosaveest=['Momentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif maxincerr==0
        filetosaveest=['EstrhoMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif Rov_sb==1
        filetosaveest=['SBMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif buildup==1
        filetosaveest=['BuildupMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif one_per_aut==0   % this is the benchmark option
        filetosaveest=['No_1_per_aut_Momentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    end
    
    % specify a markov process for individual income
    [incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,3);
    
    %just for the moment, turn this into an i.i.d. process
    
    incomevals=exp(incomevals');
    QQQ=Q1^1000;
    SS=QQQ(1,:)';
    incomevals=incomevals/(SS'*incomevals);
    cumQ=cumsum(Q1,2);
    incomevals=(1-ytaxrate)*incomevals+ytaxrate; % tax as in Krueger and Perri JET 2011
    
    
    % set minimum and maximum village size
    if coaldev==1                      
        vss_set=[1:7];
        vss_high=vss_set(end);
    elseif coaldev==0
        if buildup==0% for individual deviations, start with small village and build up, but coarser
            vss_set=[vss_high_inddev];
            lagind=1;
            vss_high=vss_set(end);% an indicator
        elseif buildup==1
            vss_set=vssbuildupvec_inddev;
            lagind=1;
            vss_high=vss_set(end);% an indicator
        end
    elseif coaldev==2;
        vss_set=vss_high_inddev;
        vss_high=vss_set(end);
    elseif coaldev==3
        vss_set=[1:10];
        vss_high=10;
    end
    
    
    % ===================================
    % iterate over discount factors
    % ===================================
    if coaldev==0
        betavec=betavec0;
    elseif coaldev==1
        betavec=betavec1;
    elseif coaldev==2
        r=0.04;
        betamax=1/(1+r)-0.015;
        betavec=linspace(0.8,betamax,16);
        Z=incomevals;
        phinat=min(Z)/(r-1);
        phivec=linspace(phinat,0,Nphi);
    elseif  coaldev==3
        betavec=betavec3;
    end
    Income0=nan(vss_high_inddev+1,T,size(betavec,2),size(rhovec,2),Nincproc);
    Csim0=nan(vss_high_inddev+1,T,size(betavec,2),size(rhovec,2),Nincproc);
    Momentsest=nan(31,N_cerr+1,size(betavec,2),1,vss_high_inddev);
    for Qbeta=1:size(betavec,2) 
        
        % ===================================
        % iterate over CRRA
        % ===================================
        for Qrho=1:size(rhovec,2)
            beta=betavec(Qbeta)
            rho=rhovec(Qrho)
            
            % limited commitment economies
            
            if coaldev==1
                lagind=[];
            else
                lagind=1;
            end
            
            P_LAG=nan(1000,1000,length(vss_set));Y_lag=nan(1000,1000,length(vss_set));
            WAUT_LAG=nan(1000,1000,length(vss_set));
            waut_lag=[]; STABLE=[]; STABLECHECK=[];
                        if coaldev==3
                Stable=nan(length(vss_set));
                stable=[]; countstable=0;
                clear Delta
            end
            for vsscount=1:length(vss_set)       % iterate over village size
                vss=vss_set(vsscount);
                vss_max=max(vss_set);
                if vsscount>1 & coaldev==1
                    lagind=lagind+(vss-(vss_set(vsscount-1)+1));
                end
                
                disp('now solving for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
                % Aiyagari economy
                if coaldev==2
                    sav_problem
                    % Full insurance with coal deviations - simple example in
                    % Genicot and Ray   
                else
                    % this is the scaling factor: if you want RoV consumption as average consumption per HH, scale rov variables
                    % by vss
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                    csim=[];  cc =[]; eeewaut =[]; eewaut =[]; ewaut =[]; waut =[]; waut_check=[]; S_check=[]; waut_check_ind=[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov =[]; P_rov =[];
                    S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
                    wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[];
                    Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
                    S_rov_aut=[]; S_rov_sb=[]; S_rov_aut_all=[];  wsb =[]; ewsb =[];
                    w =[];  phinew =[]; Yagg =[]; wint =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; bind=[];
                    
                    
                    % ===================================
                    % build the state space of ROV income
                    % ===================================
                    idio_dist=[];
                    idio_dist_aut=[];
                    idio_dist_aut_all=[];
                    
                    
                    for k1=0:vss
                        for k2=0:vss
                            for k3=0:vss
                                II=zeros(1,3);
                                if k1+k2+k3==vss
                                    % the vector of aggregate incomes from
                                    % permutations of vss households
                                    S_rov=[S_rov; k1*incomevals(1,1) + k2*incomevals(2,1) + k3*incomevals(3,1)];
                                    % the corresponding distribution of
                                    % individual income values
                                    idio_dist=[idio_dist; k1 k2 k3];
                                    
                                    % this computes the 'best' distribution of incomes
                                    % of size vss-lagind
                                    if vss>1
                                        II=zeros(1,3);
                                        temp=[k1 k2 k3];
                                        % this takes away the lowest incomes to
                                        % derive the distribution of incomes of the
                                        % 'best' coalition lagind individuals
                                        % This coalition is one smaller than the
                                        % largest sustainable coalition because it
                                        % will consist of 1 agent and m-1 guys from
                                        % the largest sustainable coalition of size
                                        % m
                                        
                                        for jj=1:lagind
                                            II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                            temp=[k1 k2 k3]-II;
                                        end
                                        idio_dist_aut=[idio_dist_aut; [k1 k2 k3]-II];
                                        S_rov_aut=[S_rov_aut; ([k1 k2 k3]-II)*incomevals];
                                        
                                        
                                        II=zeros(1,3);
                                        temp=[k1 k2 k3];
                                        % this takes away the lowest incomes to
                                        % derive the distribution of incomes of the
                                        % 'best' coalition
                                        % consisting of m individuals from only the
                                        % rest of village
                                    end
                                end
                            end
                        end
                    end
                    
                    % resort S_rov in ascending order
                    [S_rov,S_rovind]=sort(S_rov);
                    idio_dist=idio_dist(S_rovind,:);
                    
                    if vss>1
                        idio_dist_aut=idio_dist_aut(S_rovind,:);
                        S_rov_aut=S_rov_aut(S_rovind);
                    end
                    %and create S_rov_aut vectors that drops lagind individuals from
                    %the rov vector. 
                    if vss>1
                        for k=1:length(STABLE)
                            if STABLE(k)==1 %not sure, do this for all
                                for g=1:length(idio_dist)
                                    index=0;
                                    for j1=0:vss-k
                                        for j2=0:vss-k
                                            for j3=0:vss-k
                                                if j1+j2+j3==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2 & idio_dist(g,3)>=j3)
                                                    index=index+1;
                                                    temp=[j1 j2 j3];
                                                    idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp];
                                                    S_rov_aut_all(g,index,k)=[(idio_dist(g,:)-temp)*incomevals];
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                            
                        end
                    end
                    % ===================================
                    % simulate aggregate income transitions
                    % ===================================
                    tic
                    T_trans=50000;
                    aggincind=zeros(length(S_rov),length(S_rov));
                    for jjj=1:length(S_rov)
                        jjj;
                        inccount=zeros(T_trans,vss,3);
                        agginc=zeros(T_trans,3);
                        for iii=find(idio_dist(jjj,:)>0)
                            
                            for kkk=1:idio_dist(jjj,iii)
                                rng(jjj*100+iii*10+kkk,'twister');
                                shock=rand(T_trans,1);
                                rr=(find(shock<=cumQ(iii,1)));
                                inccount(rr,kkk,iii)=1;
                                rr=(find(cumQ(iii,1)<shock & cumQ(iii,2)>=shock));
                                inccount(rr,kkk,iii)=2;
                                rr=(find(cumQ(iii,2)<shock & cumQ(iii,3)>=shock));
                                inccount(rr,kkk,iii)=3;
                            end
                            
                            agginc(:,iii)=sum(incomevals(inccount(:,1:idio_dist(jjj,iii),iii)),2);
                        end
                        agginc=sum(agginc,2);
                        for jjj2=1:length(S_rov)
                            qqq=(agginc<=S_rov(jjj2)+0.0000000001).*(agginc>=S_rov(jjj2)-0.0000000001);
                            aggincind(jjj,jjj2)=sum(qqq);
                        end
                        P_rov(jjj,:)=aggincind(jjj,:)/sum(aggincind(jjj,:));
                    end
                    disp('Simulation for transitions of income distribution takes')
                    toc

                    if vss>1 && coaldev<2
                        S_rov_aut=S_rov_aut/(max(ScaleRov-lagind,1));
                        for g=vss_set(find(vss-1>=vss_set))
                            if STABLE(g)==1
                                S_rov_aut_all(:,:,g)=S_rov_aut_all(:,:,g)/(g);
                            end
                        end
                    end
                    
                    S_rov=S_rov/ScaleRov;
                    
                    % ===================================
                    % state space: combinations of indiv and village income
                    % ===================================
                    S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
                    idio_dist_rov=kron(ones(size(incomevals)),idio_dist);
                    idio_dist_aut=kron(ones(size(incomevals)),idio_dist_aut);
                    S_rov_aut=[kron(ones(size(incomevals)),S_rov_aut)];
                    if vss>1 & coaldev==1
                        for g=vss_set(find(vss-1>=vss_set))
                            if STABLE(g)==1
                                for k=1:size(S_rov_aut_all,2)
                                    S_rov_aut_all2(:,k,g)=[kron(ones(size(incomevals)),S_rov_aut_all(:,k,g))];
                                    idio_dist_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),idio_dist_aut_all(:,:,k,g))];
                                end
                            end
                        end
                        S_rov_aut_all=S_rov_aut_all2;
                        idio_dist_aut_all=idio_dist_aut_all2;
                        clear S_rov_aut_all2;
                        clear idio_dist_aut_all2;
                    end
                    
                    P=kron(Q1,P_rov);
                    
                    Y=S(:,1)+ScaleRov*S(:,2);				% possible aggregate income outcome
                    state=length(S);                        % number of states
                    agent=2;                            	% number of 'agents' in the HH vs. RoV approximation equals 2
                    STATE(vss)=state;
                    
                    
                    % ===================================
                    % Compute set of autarky values
                    % ===================================
                    % autarky values: coalitional deviations imply that the individual's
                    % outside option is given by utility of caut for 1 period, and then the
                    % expected value being in the 'best' insurance coalition of
                    % size vss-1, with most 'equal' weights, i.e.
                    % gamma=1/(vss-1) or, if this does not yield the autarky value corresponding to a coalition of vss-2, the value that is closest to it
                    for k=1:agent
                        caut(:,k)=S(:,k);
                    end
                    waut=zeros(state,agent);
                    % this is just the PDV of utility from consumption of
                    % income caut
                    eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
                    eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
                    if coaldev<2
                        % in first iteration, or for individual deviations, waut is just consumption of income
                        if vss==1 || coaldev==0
                            for j=1:state
                                waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,1);
                                % autarky values for the rest of village are always the same - perfect
                                % risk-sharing among the rest of village
                                waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                            end
                            
                        else
                            % this calculates autarky values for the coalitional
                            
                            for j=1:state
                                % indicator for the 'autarky' state in the state
                                % space of the previous iteration S_lag
                                S_lag_2=(abs(S_rov_aut(j)-S_lag(:,2))<1.0e-15) ;
                                S_lag_1=(abs(S(j,1)-S_lag(:,1))<1.0e-15);
                                S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                % for the individual, autarky value is utility of
                                % current consumption plus relevant
                                % transition
                                % probability times waut_lag, the vector of values from
                                % being in a stable, smaller insurance group with the
                                % best partners (i.e. most desirable in terms of the probability distribution this implies for the next period).
                                waut_old(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                                clear S_lag_2 S_lag_1 S_lag_ind
                                % initialize waut_temp at low number
                                waut_temp1=-100000*(ones(10000,length(STABLE))); maxindex=0;
                                waut_temp2=-100000*(ones(10000,length(STABLE)));
                                for k=1:length(STABLE)
                                    if STABLE(k)==1
                                        index=0;
                                        for j1=0:vss-k
                                            for j2=0:vss-k
                                                for j3=0:vss-k
                                                    if j1+j2+j3==vss-k && (idio_dist_rov(j,1)>=j1 & idio_dist_rov(j,2)>=j2 & idio_dist_rov(j,3)>=j3)
                                                        index=index+1     ;
                                                        S_lag_2=(abs(S_rov_aut_all(j,index,k)-S_KEEP(:,2,k))<1.0e-15);
                                                        S_lag_1=(abs(S(j,1)-S_KEEP(:,1,k))<1.0e-15);
                                                        S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                                        waut_temp1(index,k)= P_LAG(S_lag_ind,find(~isnan(P_LAG(S_lag_ind,:,k))),k)*WAUT_LAG(find(~isnan(P_LAG(S_lag_ind,:,k))),1,k);
                                                        waut_temp2(index,k)= WAUT_LAG(S_lag_ind,1,k);
                                                        clear S_lag_2 S_lag_1 S_lag_ind
                                                        if index>10000 % make sure index doesn't exceed dimension of waut_temp
                                                            stop
                                                        end
                                                    end
                                                    
                                                end
                                            end
                                        end
                                    end
                                end
                                % take the best continuation value of all
                                % sustainable coalitions - NB: this
                                % abstracts from incentives to deviate of
                                % other coalition members, we take care of
                                % that below
                                if one_per_aut==1
                                    waut(j,1)= utility(rho,caut(j,1)*(1-Pun_fac))+beta*max(max(waut_temp1));
                                    
                                elseif one_per_aut==0
                                    waut(j,1) = max(max(waut_temp2));
                                end
                                clear waut_temp1 waut_temp2
                                % for rest of village, the autarky value is just
                                % expected disc utility of current consumption
                                waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                                

                                for g=1:size(S_rov_aut_all,3)
                                    if STABLE(g)==1
                                        for k=1:size(S_rov_aut_all,2)
                                            S_lag_2=(abs(S_rov_aut_all(j,k,g)-S_LAG(:,2,g))<1.0e-15);
                                            S_lag_1=(abs(S(j,1)-S_LAG(:,1,g))<1.0e-15);
                                            S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                            clear S_lag_2 S_lag_1
                                            agentstate=[];
                                            for kk=1:length(incomevals)
                                                if idio_dist_aut_all(j,kk,k,g)>0
                                                    agentstate=[agentstate kk];
                                                end
                                            end
                                            
                                            if size(S_lag_ind)>0
                                                for kk=1:length(agentstate) %note, you only need one entry per income state, not one entry per agent!!!!!
                                                    if one_per_aut==1
                                                        waut_check_ind(j,1,k,g)=utility(rho,S(j,1))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),1,g);
                                                        waut_check(j,kk,k,g)=utility(rho,incomevals(agentstate(kk)))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),2,g);
                                                    elseif one_per_aut==0
                                                        waut_check_ind(j,1,k,g)=WAUT_LAG(S_lag_ind,1,g);
                                                        waut_check(j,kk,k,g)=WAUT_LAG(S_lag_ind,2,g);
                                                    end
                                                    
                                                    S_check(j,kk,k,g)=agentstate(kk);  
                                                    
                                                end
                                            end
                                            clear  S_lag_ind
                                        end
                                    end
                                end
                            end
                        end

                        % ====================================================
                        % Compute constrained insurance
                        % ====================================================
                        gammagrid=[lambda vss/ScaleRov*(1-lambda)];
                        gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];
                        gammagridold=gammagrid;
                        clear gammagrid
                        gammagrid=[linspace(1.2*max(S(:,1)),min(S(:,1)),m)',linspace(min(S(:,2)),1.2*max(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,1)),min(S(:,1)),m)'./linspace(min(S(:,2)),1.2*max(S(:,2)),m)'];
                        
                        % ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
                        n=length(gammagrid);
                        w=zeros(n,state,agent);
                        gammar=zeros(n,state,agent);
                        c=zeros(n,state,agent);
                        phi=zeros(n,state,agent);
                        
                        % 'first-best' consumption given aggregate resources and planner weights
                        for j=1:state
                            c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
                            c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
                        end
                        
                        cfirstbest=c;
                        %  relative planner weights
                        for j=1:state
                            for k=1:agent
                                gammar(:,j,k)=gammagrid(:,agent+k);
                            end
                        end
                        
                        % Compute the first best value functions without binding
                        % constraints - i.e. cfirstbest for ever
                        w=zeros(n,state,agent);
                        
                        dww=1;
                        wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
                        while dww>0.00001
                            w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
                            w(:,:,2)=utility(rho,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
                            
                            dww=max(max(max(abs(w-wold))));
                            wold=w;
                        end
                        wfirstbest=w;
                        
                        % continuation utility as function of planner weight and state of the
                        % economy
                        ewfirstbest(:,:,1)=wfirstbest(:,:,1)*P';
                        ewfirstbest(:,:,2)=wfirstbest(:,:,2)*P';

                        % ====================================================
                        % Start the loop with the enforcement constraints
                        % ====================================================n
                        ewsb=ewfirstbest;
                        wsb=wfirstbest;
                        options=optimset('Display','off');
                        eps=0.0000001;t=0;itgap=0;GAP=2; gap=2; gaphist=gap; gapincrease=0;
                        
                        csolold=cfirstbest;
                        if beta<0.85
                            gapcrit=gapcritset;
                            maxitgap=maxitgapset;
                        else
                            gapcrit=gapcritset*1
                            maxitgap=maxitgapset;
                        end
                        
                        while GAP >=gapcrit && itgap<=maxitgap
                            rr=[];
                            rr1=[];
                            rr2=[];
                            rrnone=[];
                            rrboth=[];    rrnoneadd=[];
                            nonconv=zeros(state,3)*nan;
                            itgap=itgap+1;
                            coalnonsust=nan(state,2);
                            count=0; countgrid=[ones(n,state)];
                            for j=1:state
                                
                                                                
                                index=[1 j];
                                rr1=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if vss>1 & coaldev==1
                                    rr1cd=[];
                                    %collect all possible coalitions
                                    for g=size(waut_check,4):-1:0 %check from the largest group down, 0 is an individual deviation
                                        if g>0 & STABLECHECK(g+1)==1
                                            
                                            for k=1:size(waut_check,3)
                                                if waut_check(j,1,k,g)~=0
                                                    coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-waut_check_ind(j,1,k,g)];
                                                    for kk=1:size(waut_check,2)
                                                        if waut_check(j,kk,k,g)~=0
                                                            coalition=[coalition utility(rho,cfirstbest(rr1,j,2))+beta*ewsb(rr1,j,2)-waut_check(j,kk,k,g)];
                                                        end
                                                    end
                                                    groupchecksolve=coalition<0;
                                                    groupchecksolve=prod(groupchecksolve,2);
                                                    groupchecksolve=find(groupchecksolve>0);
                                                    if size(groupchecksolve)>0
                                                        [j k g];
                                                    end
                                                    COALITION(:,1:size(coalition,2),k,g+1)=[coalition];
                                                end
                                                
                                            end
                                            
                                        elseif g==0
                                            COALITION(:,1,1,g+1)=utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-eeewaut(j,1);
                                        end
                                    end
                                    %check for each grid point which coalition is the
                                    %most binding (from the point of view of the
                                    %individual)
                                    for kk=1:length(rr1)
                                        for g=1:size(COALITION,4)
                                            for k=1:size(COALITION,3)
                                                rrr=find(COALITION(kk,:,k,g)~=0);
                                                if COALITION(kk,rrr,k,g)<0
                                                    value(kk,k,g)=COALITION(kk,1,k,g);
                                                    BIND(kk,k,g)=1;
                                                else
                                                    value(kk,k,g)=9999;
                                                    BIND(kk,k,g)=0;
                                                end
                                            end
                                            [CC1,II1]=min(value(kk,:,g));
                                            value2(kk,g)=value(kk,II1,g);
                                            index2(kk,g)=[II1];
                                            BIND2(kk,g)=BIND(kk,II1,g);
                                        end
                                        [CC,II]=min(value2(kk,:));
                                        
                                        gmax(kk)=II;
                                        kmax(kk)=index2(kk,II);
                                        bindmax(kk)=BIND2(kk,II);
                                        valuemax(kk)=CC;
                                        
                                    end
                                    %collect grid points that have the same binding
                                    %coalition and solve for the solution
                                    %simultaneously
                                    RRCHECK=[];
                                    for g=1:size(COALITION,4)
                                        for k=1:size(COALITION,3)
                                            rrcheck=find(kmax==k & gmax==g & bindmax==1);
                                            RRCHECK=[RRCHECK; rrcheck'];%because these grid points have been dealt with
                                            if size(rrcheck)>0
                                                
                                                bind=1;
                                                if g>1
                                                    coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-waut_check_ind(j,1,k,g-1)];
                                                    for kk=1:size(waut_check,2)
                                                        if waut_check(j,kk,k,g-1)~=0
                                                            coalitionsolve=[coalitionsolve utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)-waut_check(j,kk,k,g-1)];
                                                        end
                                                    end
                                                elseif g==1
                                                    coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-eeewaut(j,1)];
                                                end
                                                rrchecksolve=coalitionsolve<0;
                                                rrchecksolve=prod(rrchecksolve,2);
                                                rrchecksolve=find(rrchecksolve>0);
                                                [j k g];
                                                rr=rrchecksolve(end)+1;
                                                if rr>length(gammagrid)
                                                    coalnonsust(j,1)=-2;
                                                    output=[nan,nan,-50];
                                                else
                                                    clear x0
                                                    
                                                    x0(1) = cfirstbest(rr,j,1);  
                                                    if g==1
                                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                                    elseif g>1
                                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[waut_check_ind(:,1,k,g-1) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                                    end
                                                end
                                                if output(3) <=0
                                                    nonconv(j,1)=output(3);
                                                    csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                                    csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                                    gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                    wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                                    wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                                    phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                    phisol(rr1(rrcheck),j,2:agent)=0	;
                                                    
                                                else
                                                    
                                                    csol(rr1(rrcheck),j,1)=x(1);
                                                    gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                    csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                                    wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                                    wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                                    phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                    phisol(rr1(rrcheck),j,2:agent)=0	;
                                                end
                                            end
                                        end
                                        
                                    end
                                    %and this is the remainder;
                                    rr1(RRCHECK)=[];
                                    rrnoneadd=rr1;
                                    COALITION=[];
                                    coalition=[]; coalitionsolve=[]; rrchecksolve=[];
                                    clear value BIND value2 BIND2 index2 gmax kmax bindmax valuemax

                                else
                                    %in first group size iteration, i.e. for vss=1
                                    rr1cd=rr1;
                                    rr1id=[];
                                    rrnoneadd=[];
                                end

                                if size(rr1cd)>0
                                    bind=1;

                                    rr=rr1cd(end)+1;
                                    if rr>length(gammagrid)
                                        coalnonsust(j,1)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        if size(rr1cd)>0
                                            clear x0
                                            
                                            x0(1) = cfirstbest(rr1cd(1),j,1);  %redundant because we no longer need an initial guess
                                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,waut,bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                        end
                                        
                                    end
                                    if output(3) <=0
                                        nonconv(j,1)=output(3);
                                        csol(rr1cd,j,1)=ones(size(rr1cd,1),1)*caut(j,1);
                                        csol(rr1cd,j,2)=ones(size(rr1cd,1),1)*caut(j,2);
                                        gammar(rr1cd,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=ones(size(rr1cd,1),1)*waut(j,1);
                                        wsbnewsol(rr1cd,j,2)=ones(size(rr1cd,1),1)*waut(j,2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    else
                                        csol(rr1cd,j,1)=x(1);
                                        gammar(rr1cd,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                        csol(rr1cd,j,2)=(Y(j)-x(1))/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=output(1);
                                        wsbnewsol(rr1cd,j,2)=output(2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    end
                                end
                                %this becomes relevant for vss>1, these are all the left-overs that are not dealt with by joint deviations.
                                if size(rrnoneadd)>0
                                    
                                    csol(rrnoneadd,j,1:agent)=cfirstbest(rrnoneadd,j,1:agent);
                                    gammar(rrnoneadd,j,1:agent)=gammagrid(rrnoneadd,agent+1:2*agent);
                                    phisol(rrnoneadd,j,1:agent)=0;
                                    wsbnewsol(rrnoneadd,j,1)=utility(rho,cfirstbest(rrnoneadd,j,1))+beta*ewsb(rrnoneadd,j,1);
                                    wsbnewsol(rrnoneadd,j,2)=utility(rho,cfirstbest(rrnoneadd,j,2))+beta*ewsb(rrnoneadd,j,2);
                                end
                                
                                % 2. Constraint binding for agent 2 (RoV)
                                rr2=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)) ;
                                
                                if size(rr2)>0
                                    bind=2; %i.e. agent 2 has the binding constraint
                                    rr=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                                    if isempty(rr)==1
                                        coalnonsust(j,2)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        clear x0
                                        x0(1) = cfirstbest(rr(end),j,1); %redundant because we no longer need an initial guess
                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut(:,1) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);

                                    end
                                    % output(3) equals 1 if constraints
                                    % are satisfied, but -2 if not
                                    if output(3)<=0
                                        nonconv(j,2)=output(3)*100;
                                        csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                        gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                        wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                        wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                        phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                        phisol(rr2,j,1)=0;
                                    else
                                        %first check that agent 1 (in conjunction with the rest of the village)
                                        %truly does not have binding constraints
                                        %at the proposed solution
                                        if vss>1  & coaldev==1 % only need to check this once coalitional deviations become relevant
                                            RRCHECK=[];
                                            for g=size(waut_check,4):-1:1 %check from the largest group down, 0 is an individual deviation and already dealt with above
                                                if STABLECHECK(g+1)==1
                                                    %you need to add one here, because you've got the individual in this vector as well
                                                    %check the largest possible coalition first:
                                                    
                                                    for k=1:size(waut_check,3)   %this code relies on the most binding coalition being with the individual with the highest income, may no longer be true with negative serial correlation
                                                        if waut_check(j,1,k,g)~=0
                                                            coalition=[output(1)-waut_check_ind(j,1,k,g)];
                                                            for kk=1:size(waut_check,2)
                                                                if waut_check(j,kk,k,g)~=0
                                                                    coalition=[coalition output(2)-waut_check(j,kk,k,g)];
                                                                end
                                                            end
                                                            COALITION(:,1:size(coalition,2),k,g+1)=[coalition];
                                                        end
                                                    end
                                                    
                                                end
                                            end
                                            
                                            %and now check if any of these binds
                                            for g=1:size(COALITION,4)
                                                for k=1:size(COALITION,3)
                                                    rrr=find(COALITION(1,:,k,g)~=0);
                                                    if COALITION(1,rrr,k,g)<0
                                                        RRCHECK=[RRCHECK; 1];
                                                    end
                                                end
                                            end

                                            %is there any coalition for which rest of village
                                            %would need less than waut(j,2), in which case,
                                            %we know contract is not sustainable
                                            
                                            if size(RRCHECK)>0
                                                output(3)=-2;
                                            end
                                            
                                            COALITION=[];
                                            coalition=[];
                                        end
                                    end
                                    if output(3)>0
                                        
                                        csol(rr2,j,1)=x(1);
                                        gammar(rr2,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                        csol(rr2,j,2)=(Y(j)-x(1))/ScaleRov;
                                        wsbnewsol(rr2,j,1)=output(1);
                                        wsbnewsol(rr2,j,2)=output(2);
                                        phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                        phisol(rr2,j,1)=0;
                                    else
                                        nonconv(j,2)=output(3)*100;
                                        csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                        gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                        wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                        wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                        phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                        phisol(rr2,j,1)=0;
                                    end
                                end
                                
                                % 3. Constraint all slack
                                rrnone=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                
                                if size(rrnone)>0
                                    
                                    csol(rrnone,j,1:agent)=cfirstbest(rrnone,j,1:agent);
                                    gammar(rrnone,j,1:agent)=gammagrid(rrnone,agent+1:2*agent);
                                    phisol(rrnone,j,1:agent)=0;
                                    
                                    wsbnewsol(rrnone,j,1)=utility(rho,cfirstbest(rrnone,j,1))+beta*ewsb(rrnone,j,1);
                                    wsbnewsol(rrnone,j,2)=utility(rho,cfirstbest(rrnone,j,2))+beta*ewsb(rrnone,j,2);
                                end
                                
                            end

                            gapold=max(norm(wsbnewsol(:,:,1)-wsb(:,:,1)),norm(wsbnewsol(:,:,2)-wsb(:,:,2)));
                            gap1=max(max(max(abs(invutility(beta,rho,wsbnewsol(:,:,1))-invutility(beta,rho,wsb(:,:,1)))./invutility(beta,rho,wsb(:,:,1)))),max(max(abs(invutility(beta,rho,wsbnewsol(:,:,2))-invutility(beta,rho,wsb(:,:,2)))./invutility(beta,rho,wsb(:,:,2)))))/(1-beta);
                            gap2=max(max(max(abs(csol(:,:,1)-csolold(:,:,1))./csol(:,:,1))),max(max(abs(csol(:,:,2)-csolold(:,:,2))./csol(:,:,2))));
                            gap=max(gap1,gap2);
                            gaphist(itgap+1)=gap;
                            GAP=max(gaphist(max(1,length(gaphist)-convlength):end));
                            disp('gap,gapold,gap2,gap1')
                            disp([gap,gapold,gap2,gap1])
                            nonconvind(Qbeta,Qrho,vss,incproc)=mean(nonconv(find(~isnan(nonconv))));
                            Coalnonsust(Qbeta,Qrho,vss,incproc)=mean(coalnonsust(find(~isnan(coalnonsust))));
                            
                            if  (min(min(nonconv))<=0 ||min(min(coalnonsust))<0) & itgap>=3 & coaldev==1
                                itgap
                                disp('Autarky or coalition not sustainable')
                                [nonconv,S]
                                for j=1:state
                                    csol(:,j,1)=caut(j,1);
                                    csol(:,j,2)=caut(j,2);
                                end
                                % pause
                                break
                            end
                            
                            if gaphist(end)<gaphist(end-1)
                                gapincrease=1;
                            end
                            if gapincrease==0
                                csolold=csol;
                                wsb=wsbnewsol;
                            else
                                csolold=update*csol+(1-update)*csol;
                                wsb=update*wsbnewsol+(1-update)*wsbnewsol;
                                
                            end
                            ewsb(:,:,1)=(P*wsb(:,:,1)')'; 
                            ewsb(:,:,2)=(P*wsb(:,:,2)')';                           
                        end
                        
                        CRITGAP(vss)=gaphist(itgap-1);
                        
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                        
                        
                        % =============================================================
                        % calculate values of outside option for coalitional deviations
                        % =============================================================
                        % if insurance group of size vss was unsustainable, use
                        % previous value for outside option
                        
                        if coaldev==1 && (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0) && itgap>=3 && vss>1
                            % this is redundant, but shows you that if vss is
                            % unsustainable, the outside option for vss+1 remains that
                            % for vss
                            waut_lag(:,1)=waut_lag(:,1);
                            waut_lag(:,2)=waut_lag(:,2);
                            P_lag=P_lag;
                            S_lag=S_lag;
                            % if the current insurance group is not sustainable, it
                            % cannot be the outside option for the next bigger group,
                            % so increase the 'lagind' indicator by one
                            lagind=lagind+1;
                            gap=-1;
                            STABLE(vss)=0;
                            STABLECHECK=[1 STABLE];
                            % ... or if there is no previous, use just conditional exp of
                            % autarky consumption as calculated above
                        elseif (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0)  && itgap>=3 && (vss==1 || coaldev==0)
                            waut_lag(:,1)=eeewaut(:,1);
                            waut_lag(:,2)=eeewaut(:,2);
                            P_lag=P;
                            S_lag=S;
                            lagind=1;
                            gap=-1;

                            S_rov_KEEP(1:state/3,vss)=S_rov;
                            S_KEEP(1:state,:,vss)=S;
                            PHI_KEEP(:,1:state,:,vss)=phisol;
                            STABLE(vss)=1;
                            STABLECHECK=[1 STABLE];
                            IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss)=idio_dist;
                        else     
                            lagind=1;
                            % when ins group of size vss is sustainable, its value
                            % becomes outside option for next bigger group size
                            for j=1:state
                                % outside option for next larger insurance group is
                                % value with equal weights vss ...
                                eewaut(j,1)=wsb((n+1)/2,j,1);
                                eewaut(j,2)=wsb((n+1)/2,j,2);
                                
                            end
                            waut_lag=eewaut;
                            S_lag=S;
                            P_lag=P;
                            
                            S_rov_KEEP(1:state/3,vss)=S_rov;
                            S_KEEP(1:state,:,vss)=S;
                            PHI_KEEP(:,1:state,:,vss)=phisol;
                            WAUT_LAG(1:state,1:size(waut_lag,2),vss)=waut_lag;
                            P_LAG(1:state,1:state,vss)=P_lag;
                            S_LAG(1:state,1:size(S_lag,2),vss)=S_lag;

                            STABLE(vss)=1;
                            STABLECHECK=[1 STABLE]; %this includes stability for the individual as well, makes the loop above easier
                            IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss)=idio_dist;
                        end
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                    elseif coaldev==3
                        % stationary distribution of aggregate income states
                        SS_S=P^300;
                        SSS_S=SS_S(1,:);
                        SSS_Ssave(vss+1,1:length(SSS_S))=SSS_S;
                        sizeSSS_Ssave(vss+1)=length(SSS_S);
                        % value from consuming average income
                        if vss==1
                            sizeSSS_Ssave(1)=3;
                            SSS_Ssave(1,1:3)=SS;
                            V_FI(1,1:3)=(inv(diag(ones(length(incomevals),1))-beta*Q1)*utility(rho,incomevals));
                            P_LAG(1:3,1:3,vss)=Q1;
                            Y_lag(1,1:3,1)=incomevals';
                        end
                        V_FI(vss+1,1:length(S))=(inv(diag(ones(length(S),1))-beta*P)*utility(rho,Y/(vss+1)))';
                        
                        % vvss are the number of high income individuals
                        Y_lag(vss+1,1:length(Y))=Y;
                        P_LAG(1:length(S),1:length(S),vss+1)=P;
                        % vvss are the number of high income individuals,
                        % required to correspond to stable group sizes. NB:
                        % individuals and pairs are always stable
                        for vvss=[1,stable'+1]
                            ind1=find(abs(Y-(vvss)*max(incomevals)-(vss+1-vvss)*min(incomevals))<0.00000001); % this is the indicator of the state today with vvss high and vss+1-vvss low income individuals - can be several
                            ind2=find(abs(Y_lag(vvss,:)-(vvss)*max(incomevals))<0.000000001); % this is the indicator of the State in a group of vvss individuals associated with only high income individuals
                            % this is the utility loss from making transfers today when
                            % vvss out of vss individuals have high income and vss-vvss
                            % have minimum income
                            deltav(vss,vvss)=(utility(rho,max(incomevals))-utility(rho,(vvss*max(incomevals)+(vss+1-vvss)*min(incomevals))/(vss+1)));
                           
                            % the difference in continuation values from having full
                            % insurance with vss vs vvss individuals
                            % NB: use ind1(2), which corresponds to the
                            % individuual having high income
                            check1(vss,vvss)=beta*(P(ind1(end),:)*V_FI(vss+1,1:sizeSSS_Ssave(vss+1))');
                            check2(vss,vvss)=beta*(P(ind1(1),:)*V_FI(vss+1,1:sizeSSS_Ssave(vss+1))');
                            deltaV(vss,vvss)=beta*(P(ind1(end),:)*V_FI(vss+1,1:sizeSSS_Ssave(vss+1))'-P_LAG(ind2(1),1:sizeSSS_Ssave(vvss),vvss)*V_FI(vvss,1:sizeSSS_Ssave(vvss))');
                            % if this is negative for vvss, and the group with vvss is stable, then full insurance with vss individuals
                            % is sustainable against a deviation with vvss individuals
                            Delta(vss,vvss)=deltav(vss,vvss)-deltaV(vss,vvss)
                            
                        end
                        Stable(vss,stable)=(Delta(vss,stable)<0);
                        if sum(Delta(vss,:)<0)==length(stable)+1
                            countstable=countstable+1;
                            stable(countstable,1)=vss;
                        end
                        clear yind  y yy dy aggyy aggy
                    end
                end
            end
            
            clear P_LAG WAUT_LAG
            %==========================================================
            % Simulate the economy: CD only for the largest group size;
            % ID for all
            % =========================================================
            % This lists all possible combinations of individual income
            % values (IDIO_DIST), and the aggregate states (pair of individual income and ROV-income) associated
            % with every individual income value
            
            if coaldev==1                       % for coalitional deviations, all the stable ones
                stable=(find(STABLE==1));
                vss_set_sim=stable;
            elseif coaldev==0                                % for individual deviations, buildup (but with holes)
                vss_set_sim=vss_set;
            elseif coaldev==2
                vss_set_sim=vss_set;
            elseif coaldev==3
                vss_set_sim=stable';
            end
            
            for vss=[vss_set_sim]
                if coaldev<2
                    state=STATE(vss);
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                    
                    S=S_KEEP(1:state,:,vss);
                    S_rov=S_rov_KEEP(1:state/3,:,vss);
                    idio_dist=IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss);
                    phisol=PHI_KEEP(:,1:state,:,vss);
                    IDIO_DIST=[ones(state/3,1)*[1 0 0] + idio_dist; ones(state/3,1)*[0 1 0] + idio_dist; ones(state/3,1)*[0 0 1]+idio_dist];
                    
                    distribution=size(IDIO_DIST,1);
                end
                
                Nvill=ceil((vss_high_inddev+1)/(vss+1))
                var_dyagglog=zeros(Nvill*(vss+1),vss_high);var_dyvillagglog=zeros(Nvill*(vss+1),vss_high);vardaggylog=zeros(Nvill*(vss+1),vss_high);betayc=zeros(2,Nvill*(vss+1));vardcshare=zeros(Nvill*(vss+1),vss_high);vardcshare_dylow=zeros(Nvill*(vss+1),vss_high);vardcshare_dyhigh=zeros(Nvill*(vss+1),vss_high);varclog=zeros(Nvill*(vss+1),vss_high);vardclog=zeros(Nvill*(vss+1),vss_high);varylog=zeros(Nvill*(vss+1),vss_high);vardylog=zeros(Nvill*(vss+1),vss_high);varclog_cond=zeros(Nvill*(vss+1),vss_high);varylog_cond=zeros(Nvill*(vss+1),vss_high);
                varc_yhigh_cond=zeros(Nvill*(vss+1),vss_high);varc_ylow_cond=zeros(Nvill*(vss+1),vss_high);vardclog_cond=zeros(Nvill*(vss+1),vss_high);vardylog_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high);
                vardc_dypos_cond=zeros(Nvill*(vss+1),vss_high);
                vardc_dyneg_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high); varc_ylow=zeros(Nvill*(vss+1),vss_high);vardc_dypos=zeros(Nvill*(vss+1),vss_high); vardc_dyneg=zeros(Nvill*(vss+1),vss_high);
                
                betadcdy_agg=zeros(3,Nvill*(vss+1));   betadcdy=zeros(2,vss_high);  betadcdy_high=zeros(1,Nvill*(vss+1));betadcdy_low=zeros(1,Nvill*(vss+1));relcvar=zeros(Nvill*(vss+1),vss_high); reldcvar=zeros(Nvill*(vss+1),vss_high); reldc0var=zeros(Nvill*(vss+1),vss_high); relcvar_cond=zeros(Nvill*(vss+1),vss_high); reldcvar_cond=zeros(Nvill*(vss+1),vss_high);
                beta_cond=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_cond=zeros(4,Nvill*(vss+1),vss_high); beta_yhigh_cond=zeros(4,Nvill*(vss+1),vss_high) ;
                beta_agg=zeros(2,Nvill*(vss+1),vss_high); beta_yhigh_agg=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_agg=zeros(2,Nvill*(vss+1),vss_high);

                %=============================================
                % Start simulation
                % ============================================
                % number of insurance groups equals village size
                % vss_high+1 divided by size of ins group vss+1
                
                stream = RandStream.getGlobalStream;
                Npanel=Nvill*(vss+1);
                csim0=nan(Nvill,vss+1,T);
                income0=nan(Nvill,vss+1,T);
                income0ind=nan(Nvill,vss+1,T);
                phisim0=nan(Nvill,vss+1,T);
                Yvill=nan(1,T);
                Yagg=nan(Nvill,T);
                Ssim=nan(Nvill,vss+1,T);
                gammarfun=nan(Nvill,vss+1,T);
                phinew=nan(Nvill,vss+1,T);
                indSjj=nan(T,vss+1);
                
                disp('simulating')
                
                disp('Reminder: now simulating for village,coaldev,incproc')
                disp([village coaldev incproc])
                disp('now simulating for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
                
                tic
                
                %initial state
                t=1;
                % initial values of income and consumption (set =
                % income)
                rng(3,'twister');
                shockinittemp=randi(3,Nvill*(vss+1),1);
                shockinit=reshape(shockinittemp,Nvill,vss+1);
                income0ind(:,:,1)=shockinit;
                income0(:,:,1)=incomevals(income0ind(:,:,1));
                csim0(:,:,1)=income0(:,:,1);
                Yagg(:,t)=sum(income0(:,:,t),2);
                IDIO_DIST_sim=[];
                
                while t<T
                    t=t+1;
                    %simulate next period income taking into account
                    % persistence across time
                    rng(t,'twister');
                    shocktemp=rand(Nvill*(vss+1),1);
                    shock=reshape(shocktemp,Nvill,vss+1);
                    for k=1:vss+1
                        rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                        income0ind(rr,k,t)=1;
                        rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                        income0ind(rr,k,t)=2;
                        rr=(find(cumQ(income0ind(:,k,t-1),2)<shock(:,k) & cumQ(income0ind(:,k,t-1),3)>=shock(:,k)));
                        income0ind(rr,k,t)=3;
                        income0(:,k,t)=incomevals(income0ind(:,k,t))  ;
                    end
                    
                    Yagg(:,t)=sum(income0(:,:,t),2);
                end
                
                if coaldev<2
                    
                    t=1;
                    while t<T
                        t=t+1;
                        gammarfun(:,:,t)=(Yagg(:,t-1)*ones(1,vss+1)-csim0(:,:,t-1))./csim0(:,:,t-1)/ScaleRov;
                        gammarfun(:,:,t)=min(gammarfun(:,:,t),max(gammagrid(:,4))*ones(Nvill,vss+1));

                        for ivill=1:Nvill
                            IDIO_DIST_sim=find(((IDIO_DIST(:,1,1)==(sum(income0ind(ivill,:,t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(ivill,:,t)==2))).*(IDIO_DIST(:,3,1)==(sum(income0ind(ivill,:,t)==3))))==1)';
                            for k=1:vss+1
                                % choose for individual k the aggregate state
                                % that applies
                                 indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)==income0(ivill,k,t));
                                % check that aggregate and individual income0s are
                                % consistent with the states caculated above
                                if abs(Yagg(ivill,t)-(S(indSjj(t,k),1)+ScaleRov*S(indSjj(t,k),2)))>0.000001
                                    stop
                                end
                            end
                            Ssim(ivill,:,t)=indSjj(t,:);
                        end
                        for k=1:vss+1
                            BLA=interp1(gammagrid(:,4),phisol(:,Ssim(:,k,t),1),gammarfun(:,k,t));
                            phinew(:,k,t)=diag(BLA);
                        end
                        
                        gammaraux=((1+phinew(:,:,t))./(1+phinew(:,1,t)*ones(1,vss+1))).^(1/rho).*csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,vss+1));
                        csim0(:,:,t)=(gammaraux./(sum(gammaraux,2)*ones(1,vss+1))).*(Yagg(:,t)*ones(1,vss+1));                        
                    end
                    AA=[];
                    BB=[];
                    CC=[];
                    for t=1:T
                        AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    csim0=AA;
                    income0ind=CC;
                    income0=BB;
                elseif coaldev==2
                    
                    AA=[];
                    BB=[];
                    CC=[];
                    
                    for t=1:T
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    
                    income0ind=CC;
                    income0=BB;
                    sim_save_prob
                elseif coaldev==3
                    t=0;
                    while t<T
                        t=t+1;
                        
                        csim0(:,:,t)=(Yagg(:,t)*ones(1,vss+1))./(vss+1);
                    end
                    
                    AA=[];
                    BB=[];
                    CC=[];
                    for t=1:T
                        AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    csim0=AA;
                    income0ind=CC;
                    income0=BB;
                end
                disp('simulating takes')
                tsim(Qbeta,Qrho)=toc
                tic
                disp('reshaping takes')
                toc

                % =========================================================
                % calculate moments of consumption and income growth
                % =========================================================
                csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
                rng(3, 'twister'); cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
                rng(4, 'twister'); incerr=randn(size(csim0));   Momentsest1=nan(30,N_cerr); ind=[];
                % loop over size of cons measurement error
                for kkk=1:size(csim0,1)
                    tempvar(kkk)=var(log(csim0(kkk,:)));
                end
                Tempvar=mean(tempvar);
                for ic=1:N_cerr+1
                    csim=exp(log(csim0)+((ic-1)/N_cerr*maxcerr/(1-(ic-1)/N_cerr*maxcerr)*Tempvar)^0.5*cerr);
                    logcsim=log(csim);
                    logcsim0=log(csim0);
                    logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
                    logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
                    % add inc meas error back in
                    income=exp(log(income0)+varnu(incproc)^0.5*incerr);
                    logincome=log(income);
                    logincome0=log(income0);
                    logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
                    income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
                    
                    for ivill=1:Nvill
                        aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logaggvillincome=log(aggvillincome);
                        logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    end
                    aggincome=sum(income,1);
                    
                    % === discard first Tinit-1 periods
                    % this chooses log or level consumption for
                    % moments
                    if momentclevinlog==1
                        ccc=logcsim(:,Tinit:T);
                        ccc_cond=logcsim_meandev(:,Tinit:T);
                        ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
                    else
                        ccc=(csim(:,Tinit:T));
                        ccc_cond=(csim_meandev(:,Tinit:T)); %NB this doesn't actually exist...
                    end
                    yyy=logincome(:,Tinit:T);
                    yyy_cond=logincome_meandev(:,Tinit:T);
                    yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
                    
                    cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
                    cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
                    cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
                    
                    yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
                    yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
                    aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
                    vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
                    
                    % calculate moments for each individiual
                    for ind=1:size(csim,1)
                        dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
                        dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
                        dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
                        dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
                        
                        % conditioning on income increases, and income
                        % decreases
                        % NB: we allocate no change to decrease!!
                        rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
                        rrhc=find(yy(ind,:)>0);
                        rrlc=find(yy(ind,:)<=0);
                        rrhc_cond=find(yy_cond(ind,:)>0);
                        rrlc_cond=find(yy_cond(ind,:)<=0);
                        rrhc_groupcond=find(yy_groupcond(ind,:)>0);
                        rrlc_groupcond=find(yy_groupcond(ind,:)<=0);
                        
                        %THE VARIANCES
                        vardclog(ind,vss)=var(cc(ind,:));
                        vardylog(ind,vss)=var(yy(ind,:));
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %conditional on village income!
                        vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
                        vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
                        vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
                        vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));
                        
                        %conditional on group income!
                        %Step 1: create your residuals
                        vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
                        vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
                        vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
                        vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %and the relative variances
                        reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));
                        reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
                        reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);
                        
                        
                        %THE BETA COEFFICIENTS
                        
                        %1) Unconditional
                        %risk-sharing whole sample
                        covtemp=cov(cc(ind,:),aggyy(1,:));
                        betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
                        covtemp=cov(cc(ind,:),yy(ind,:));
                        betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
                        [betadcdytemp,temp1,rtemp,temp2,statstemp]=regress(cc(ind,:)',[ones(1,size(yy,2))',aggyy(1,:)',yy(ind,:)']);
                        rrstat(ind,vss)=rtemp(2);
                        betadcdyalt(ind,vss)=betadcdytemp(2);
                        RR2_dcdy(ind,vss)=statstemp(1);
                        clear rtemp betadcdytemp statstemp
                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY=cc(ind,:)';
                        yy_hc(ind,:)=zeros(1,size(cc,2));
                        yy_hc(ind,rrhc)=yy(ind,rrhc);
                        oneone=ones(1,size(cc,2));
                        ones_hc(ind,:)=zeros(1,size(cc,2));
                        ones_hc(ind,rrhc)=oneone(rrhc);
                        XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
                        beta_yhigh(:,ind,vss)=regress(YY,XX);
                        betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
                        betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
                        
                        %2) conditional on village income
                        covtemp=[];
                        X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
                        YY=cc_cond(ind,:)';
                        beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        
                        X=[ones(size(aggyy(1,:),2),1),aggyy'];%
                        YY=cc(ind,:)';
                        beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        
                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY_cond=cc_cond(ind,:)';
                        yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
                        oneone=ones(1,size(cc_cond,2));
                        ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
                        XX_cond=[ones(1,size(cc_cond,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
                        beta_yhigh_cond(:,ind,vss)=regress(YY_cond,XX_cond);

                        %3) conditional on group income
                        covtemp=[];
                        X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
                        YY=cc_groupcond(ind,:)';
                        beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        
                        X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
                        YY=cc(ind,:)';
                        beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        
                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY_groupcond=cc_groupcond(ind,:)';
                        yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
                        oneone=ones(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
                        XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
                        beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);
                    end
                    var_dyagglog(1,vss)=var(aggyy);
                    R2_dcdy=mean(RR2_dcdy(1:size(csim,1),vss));
                    % vector of moments
                    % col 1-5 are parameters
                    Momentsest1(1:5,ic)=[vss,beta,rho,rho1vec(incproc),varnu(incproc)]';
                    % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
                    % (8:9), conditonal on group income (10:11)
                    Momentsest1(6:11,ic)=[mean(vardclog(1:size(csim,1),vss))/mean(vardylog(1:size(csim,1),vss)), mean(reldcvar(1:size(csim,1),vss))...
                        mean(vardclog_cond(1:size(csim,1),vss))/mean(vardylog_cond(1:size(csim,1),vss)), mean(reldcvar_cond(1:size(csim,1),vss))...
                        mean(vardclog_groupcond(1:size(csim,1),vss))/mean(vardylog_groupcond(1:size(csim,1),vss)), mean(reldcvar_groupcond(1:size(csim,1),vss))]';
                    % 12-14: beta, beta high, beta low unconditional
                    Momentsest1(12:14,ic)=[mean(betadcdy(1:size(csim,1),vss)),mean(betadcdy_high(1:size(csim,1),vss)),mean(betadcdy_low(1:size(csim,1),vss))]';
                    % 13_16: beta, beta high, beta low conditional on
                    % village income
                    Momentsest1(15:17,ic)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(4,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)]';
                    % 18:20: beta, beta high, beta low conditional on
                    % group income
                    Momentsest1(18:20,ic)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(4,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)]';
                    %21:22 are variances of dclog, dylog
                    Momentsest1(21:22,ic)=[mean(vardclog(1:size(csim,1),vss)),mean(vardylog(1:size(csim,1),vss))];
                    %23:28, relative variance dc/dy (dy>0) and var dc/dy(dy<0);
                    %unconditinoal 23:24, conditional on village income 25:26,
                    %conditional on group income 27:28
                    Momentsest1(23:24,ic)=[mean(vardc_dypos(1:size(csim,1),vss));mean( vardc_dyneg(1:size(csim,1),vss))];
                    Momentsest1(25:26,ic)=[mean(vardc_dypos_cond(1:size(csim,1),vss));mean(  vardc_dyneg_cond(1:size(csim,1),vss))];
                    Momentsest1(27:28,ic)=[mean(vardc_dypos_groupcond(1:size(csim,1),vss));mean( vardc_dyneg_groupcond(1:size(csim,1),vss))];
                    Momentsest1(30,ic)=[mean( vardylog_cond(1:size(csim,1),vss))];
                    Momentsest1(31,ic)=mean(RR2_dcdy(1:size(csim,1),vss));
                    %and the critical gap
                    if coaldev<2
                        Momentsest1(29,ic)=CRITGAP(vss);
                    elseif coaldev==2
                        Momentsest1(29,ic)=phi;
                    elseif coaldev==3
                        Momentsest1(29,ic)=nan;
                    end
                end
                
                Momentsest(:,:,Qbeta,Qrho,vss)=[Momentsest1];
                %  Income0(:,:,Qbeta,Qrho,vss)=income0(1:vss_high_inddev+1,:);
                %  Csim0(:,:,Qbeta,Qrho,vss)=csim0(1:vss_high_inddev+1,:);
                toc
            end
            
            % pause
            if buildup==0
                cd(outputpath)
                %save(filetosaveest,'Momentsest');
                cd(programpath)
            elseif buildup==1
                cd(outputpath)
                %save(filetosaveest,'Momentsest');
                cd(programpath)
            end
        end
    end
    
    cd(outputpath)
    if  buildup==0 && Rov_sb==0
        %eval(['save ', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(maxincerr==0) '.mat'])
    elseif Rov_sb==1
        %eval(['save SB', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(maxincerr==0) '.mat'])
    elseif buildup==1
        %eval(['save Buildup', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(maxincerr==0) '.mat'])
    end
    cd(programpath)
end